#include "../Species/Species1D.h"
#include "../PhysicalField/PhysicalField1D.h"
#include "../Field/Field1D.h"
#include "CellParticles.h"

Species1D::Species1D(PicParams* params_, ParallelMPI* pmpi_, int ispec) :
    Species(params_, pmpi_, ispec),
    //because PartBoundCond<dim> has no default constructor
    part_bound_cond(params_, pmpi_, ispec)
{
    init_particles();
}

Species1D::~Species1D()
{

}


void Species1D::init_particles()
{
    if(species_param.species_type == "fluid")
    {
        return;
    }
    
    particles_.resize(params->dims_space[0]);

    vector<double> cell_min(1);

    for(int i = 0; i < params->dims_space[0]; i++)
    {
        particles_[i] = new CellParticles<1>();
        cell_min[0] = params->local_min[0] + i * params->cell_length[0];
        particles_[i]->initialize(params, species_param, cell_min, params->cell_length);
    }
}


void Species1D::move(PhysicalField* fields)
{
    //cout<<"move particle 1111 "<<endl;
    LocalFields E_cell_fields[2], B_cell_fields[2];
    Field1D<double>* Ex1D = static_cast<Field1D<double>*>(fields->Ex_);
    Field1D<double>* Bx1D = static_cast<Field1D<double>*>(fields->Bx_);

    B_cell_fields[0].x = params->externB[0];
    B_cell_fields[1].x = params->externB[0];
    B_cell_fields[0].y = params->externB[1];
    B_cell_fields[1].y = params->externB[1];
    B_cell_fields[0].z = params->externB[2];
    B_cell_fields[1].z = params->externB[2];

    E_cell_fields[0].y = 0.0;
    E_cell_fields[1].y = 0.0;
    E_cell_fields[0].z = 0.0;
    E_cell_fields[1].z = 0.0;

    for(int i = 0; i < params->dims_space[0]; i++)
    {
        int ii = i + params->oversize[0];
        E_cell_fields[0].x = (*Ex1D)(ii);
        E_cell_fields[1].x = (*Ex1D)(ii+1);

        if(species_param.charge == 0.0)
        {
            push_neutral(particles_[i], E_cell_fields, B_cell_fields);
        }
        else
        {
            push(particles_[i], E_cell_fields, B_cell_fields);
        }

    }


    //sort_particles();

    //correct_particles();

}

void Species1D::project(PhysicalField* fields)
{
    //cout<<"project begen ========="<<endl;
    Field1D<double>* rho1D = static_cast<Field1D<double>*>(fields->rho_s[species_number]);

    //correct_particles();

    for(int i = 0; i < params->dims_space[0]; i++)
    {
            double rhoQ_cell_fields[2];

            rhoQ_cell_fields[0] = 0.0;
            rhoQ_cell_fields[1] = 0.0;

            //if(particles_[i]->size() == 0) cout<<"particle number = 0 "<<i<<endl;
            for(auto &particle : particles_[i]->data_)
            {
                double xpn = particle.position[0] * params->cell_length_inv[0];

                int ip_ = floor(xpn);

                //Declaration and calculation of the coefficient for interpolation
                double delta;
                double coeffxp[2];

                delta   = xpn - (double)ip_;
                coeffxp[0] = 1.0 - delta;
                coeffxp[1] = delta;

                for (int iloc=0 ; iloc<2 ; iloc++)
                {
                    rhoQ_cell_fields[iloc] += coeffxp[iloc];
                }
            }
            int ii = i + params->oversize[0];
            (*rho1D)(ii)      += rhoQ_cell_fields[0];
            (*rho1D)(ii+1)    += rhoQ_cell_fields[1];
    }

    //cout<<"project begen1111 ========="<<endl;
    for(int i = 0; i < params->dims_field[0]; i++)
    {
        (*rho1D)(i) *= species_param.weight;
        //if((*rho2D)(i,j) < 0.8e-18) cout<<i<<"  "<<j<<"  "<< (*rho2D)(i,j)<<endl;
    }

    //int imax = params->dims_field[0] - 1;
    //int jmax = params->dims_field[1] - 1;
    //cout<<(*rho2D)(0,0)<<"  "<<(*rho2D)(0,1)<<"  "<<(*rho2D)(1,1)<<"  "<<(*rho2D)(1,2)<<"  "<<(*rho2D)(2,2)<<"  "<<endl;
    //cout<<(*rho2D)(imax,jmax)<<"  "<<(*rho2D)(imax,jmax - 1)<<"  "<<(*rho2D)(imax - 1,jmax - 1)<<"  "<<(*rho2D)(imax - 1,jmax - 2)<<"  "<<(*rho2D)(imax - 2,jmax - 2)<<"  "<<endl;
    //cout<<"project begen22222 ========="<<endl;
}


void Species1D::push(CellParticles<1>* cell_particles, const LocalFields E_cell_fields[2], const LocalFields B_cell_fields[2])
{
    //double charge_over_mass = static_cast<double>(particles.charge(ipart))*one_over_mass_;
    double um0, um1, um2, up0, up1, up2, p0dot, p1dot, p2dot;
    double alpha, inv_det_t, t0, t1, t2, t0_square, t1_square, t2_square, s0, s1, s2;
    double t0t1, t1t2, t2t0;
    double p0sm, p1sm, p2sm;
    double dts2, charge_over_mass;
    LocalFields E_local, B_local;

    dts2 = params->dt * 0.5;
    charge_over_mass = species_param.charge / species_param.mass;
    alpha = charge_over_mass * dts2;

    for(auto &particle : cell_particles->data_)
    {
        interpolate(E_cell_fields, E_local, particle.position);
        interpolate(B_cell_fields, B_local, particle.position);

        um0 = particle.velocity[0] + alpha * E_local.x;
        um1 = particle.velocity[1] + alpha * E_local.y;
        um2 = particle.velocity[2] + alpha * E_local.z;

        //Rotation in the magnetic field
        t0    = alpha * B_local.x;
        t1    = alpha * B_local.y;
        t2    = alpha * B_local.z;
        t0_square   = t0*t0;
        t1_square   = t1*t1;
        t2_square   = t2*t2;
        inv_det_t = 2.0/(1.0+t0_square+t1_square+t2_square);
        s0    = t0 * inv_det_t;
        s1    = t1 * inv_det_t;
        s2    = t2 * inv_det_t;

        p0dot = um0 + um1 * t2 - um2 * t1;
        p1dot = um1 + um2 * t0 - um0 * t2;
        p2dot = um2 + um0 * t1 - um1 * t0;

        up0 = um0 + p1dot * s2 - p2dot * s1;
        up1 = um1 + p2dot * s0 - p0dot * s2;
        up2 = um2 + p0dot * s1 - p1dot * s0;

        //half-acceleration in the electric field
        p0sm = up0 + charge_over_mass * E_local.x * dts2;
        p1sm = up1 + charge_over_mass * E_local.y * dts2;
        p2sm = up2 + charge_over_mass * E_local.z * dts2;

        particle.velocity[0] = p0sm;
        particle.velocity[1] = p1sm;
        particle.velocity[2] = p2sm;

        //move the particle
        for(int i = 0 ; i < params->n_dim_field ; i++)
        {
            //particle.position_old[i]  = particle.position[i];
            particle.position[i] += params->dt * particle.velocity[i];
        }
    }

}


void Species1D::push_neutral(CellParticles<1>* cell_particles, const LocalFields E_cell_fields[2], const LocalFields B_cell_fields[2])
{
    //log_info<<"push_neutral test "<<species_param.species_name;
    for(auto &particle : cell_particles->data_)
    {
        for(int i = 0 ; i < params->n_dim_field ; i++)
        {
            //particle.position_old[i]  = particle.position[i];
            particle.position[i] += params->dt * particle.velocity[i];
        }
    }

}


void Species1D::sort_particles()
{
    int new_index[1];

    psi_particles_.clear();
    exchange_particles_.clear();

    for(int i = 0; i < params->dims_space[0]; i++)
    {
        auto cell_particles = particles_[i];
        for(int i_particle = 0; i_particle < cell_particles->size(); i_particle++)
        {
            //keep_particle = :  0 delete; 1 MPI change; 2 move to other cell; 3 still in the cell
            auto &particle = cell_particles->data_[i_particle];
            new_index[0] = floor( (particle.position[0] - params->local_min[0]) * params->cell_length_inv[0] );

            int keep_particle = -1;
            if(new_index[0] == i)
            {
                keep_particle = 3;
            }
            else if(new_index[0] >= 0 && new_index[0] < params->dims_space[0])
            {
                keep_particle = 2;
            }
            else
            {
                double old_position = particle.position[0];
                keep_particle = part_bound_cond.apply(particle);

                if(keep_particle == 3)
                {
                    new_index[0] = floor( (particle.position[0] - params->local_min[0]) * params->cell_length_inv[0] );
                    if(new_index[0] == i)
                    {
                        keep_particle = 3;
                    }
                    else if(new_index[0] >= 0 && new_index[0] < params->dims_space[0])
                    {
                        keep_particle = 2;
                    }
                    else
                    {
                        cout<<"particle boudnary error  "<<new_index[0]<<"  "<<particle.position[0]<<"  "<<old_position<<endl;
                    }
                }

            }

            if (keep_particle == 0)
            {
                psi_particles_.push_back(particle);
                cell_particles->add_particle_delete_list(i_particle);
            }
            else if(keep_particle == 1)
            {
                exchange_particles_.push_back(particle);
                cell_particles->add_particle_delete_list(i_particle);
            }
            else if(keep_particle == 2)
            {
                particles_[new_index[0]]->add_new_particle(particle);
                cell_particles->add_particle_delete_list(i_particle);
            }
            else if(keep_particle == -1)
            {
                log_error<<"error: keep_particle = -1";
            }

            //to be tested: how about performing psi_particles_.push_back(particle) and exchange_particles_.push_back(particle)
            //              in another loop
        }
    }
}


void Species1D::correct_particles()
{
    for(int i = 0; i < params->dims_space[0]; i++)
    {
        auto cell_particles = particles_[i];
        cell_particles->delete_useless_particles();
    }
}


